function [ Price, Error ] = q4cv2(  )
%Q4 Summary of this function goes here
%   Implementation of MC method using Standard constructions
%   Apply Control Variate Method using European Geometric average Call  

% Parameters
S0 = 100; K = 100; r = 0.05; sigma = 0.3;
M = 10000; N = 32; m = 100; % in sample simulations
T = 1;
dt = 1/N;
disc = exp(-r*T);

% Standard Construction
Wt = randn(M + m,N);
factor = exp((r-0.5*sigma^2)*dt + sigma*sqrt(dt)*Wt);
St = S0*cumprod(factor,2);

% Now in sample estimation
d_Payoff = disc*max(mean(St(1:M,:),2)-K,0);
% CV is the European geometric average call option
CV = disc*max(prod(St(1:M,:),2).^dt-K,0);

% Now calculate the analytical mean
miu_hat = 0.5*(r - sigma^2/2)*(1+1/N);
sig_hat2 = sigma^2/3*(1+1/N)*(1+0.5/N);
S_hat = S0*exp(sig_hat2/2 + miu_hat - r);
[call, put] =  blsprice(S_hat,K,r,T,sqrt(sig_hat2));

% Get b_hat as best linear fit slope
b = polyfit(CV,d_Payoff,1); 
bhat = b(1);

% Now out of sample pricing
d_Payoff = disc*max(mean(St(M+1:M+m,:),2)-K,0);
CV = disc*max(prod(St(M+1:M+m,:),2).^dt-K,0);
P_CV = d_Payoff - bhat*(CV-call);

Price = mean(P_CV);
Error = sqrt(var(P_CV));


end

